
#=====================================================
## Title:        Is Clientelism (Only) for the Poor?
## Authors:      Cammett, Melani, Dominika Kruszewska-Eduardo, Christiana Parreira, & Sami Atallah
## Description:  data analysis
## Date:         February 8, 2024
#=====================================================

rm(list=ls())

ipak <- function(pkg){  new.pkg <- pkg[!(pkg %in% installed.packages()[, "Package"])]
if(length(new.pkg)) install.packages(new.pkg, dependencies = TRUE)
sapply(pkg, require, character.only = TRUE)}

packages <- c("memisc", "cregg", "car", "ggplot2", "reshape2", "stargazer", "multiwayvcov", "tikzDevice",
              "lmtest", "plm", "dplyr", "tibble", "tidyr", "coefplot", "xtable") 
ipak(packages)
options(stringsAsFactors = FALSE)

survey <- read.csv("data_lebanon.csv")

#==================== data cleaning and reshaping for conjoint analysis

survey$respondent_sect <- recode(survey$confession, "24" = "None", "25"="Christian", "26"="Sunni",
                                 "27"="Shia", "28"="Druze", "29"="Christian", "30"="Christian", "31"="Christian",
                                 "32"="Christian", "33"="Christian", "34"="Christian", "35"="Christian",
                                 "36"="Alawite", "37"="Other")

survey$low_ses <- ifelse(survey$household_income==7 | survey$household_income==8, 1, 0)
survey$high_ses <- ifelse(survey$household_income==5 | survey$household_income==6, 1, 0)

year <- c(1952:1996)
survey$dob <- NULL
for (i in 1:length(as.numeric(survey$age))){
  survey$dob[i] <- year[as.numeric(survey$age)[i]]}

survey$years_old <- 2017-survey$dob

survey$respondent_sex <- ifelse(survey$sex==1, 1, 0)

cols1 = survey[,c("rally_1_1", "vote_choice", grep("Rd_1_A_",names(survey),value=TRUE),"V1")]
cols2 = survey[,c("rally_1_2", "vote_choice", grep("Rd_1_B_",names(survey),value=TRUE),"V1")]
cols3 = survey[,c("rally_2_1","vote_choice2", grep("Rd_2_A_",names(survey),value=TRUE),"V1")]
cols4 = survey[,c("rally_2_2","vote_choice2", grep("Rd_2_B_",names(survey),value=TRUE),"V1")]
cols5 = survey[,c("rally_3_1","vote_choice3", grep("Rd_3_A_",names(survey),value=TRUE),"V1")]
cols6 = survey[,c("rally_3_2","vote_choice3", grep("Rd_3_B_",names(survey),value=TRUE),"V1")]
cols7 = survey[,c("rally_4_1", "vote_choice4", grep("Rd_4_A_",names(survey),value=TRUE),"V1")]
cols8 = survey[,c("rally_4_2","vote_choice4", grep("Rd_4_B_",names(survey),value=TRUE),"V1")]

cols1$selected = ifelse(cols1$vote_choice==1,1,0)
cols2$selected = ifelse(cols2$vote_choice==2,1,0)
cols3$selected = ifelse(cols3$vote_choice2==1,1,0)
cols4$selected = ifelse(cols4$vote_choice2==2,1,0)
cols5$selected = ifelse(cols5$vote_choice3==1,1,0)
cols6$selected = ifelse(cols6$vote_choice3==2,1,0)
cols7$selected = ifelse(cols7$vote_choice4==1,1,0)
cols8$selected = ifelse(cols8$vote_choice4==2,1,0)

cols1$selected_r = ifelse(cols1$rally_1_1==1,1,0)
cols2$selected_r = ifelse(cols2$rally_1_2==2,1,0)
cols3$selected_r = ifelse(cols3$rally_2_1==1,1,0)
cols4$selected_r = ifelse(cols4$rally_2_2==2,1,0)
cols5$selected_r = ifelse(cols5$rally_3_1==1,1,0)
cols6$selected_r = ifelse(cols6$rally_3_2==2,1,0)
cols7$selected_r = ifelse(cols7$rally_4_1==1,1,0)
cols8$selected_r = ifelse(cols8$rally_4_2==2,1,0)

names(cols1) = c("rally", "vote", "religion","clientelism","security","experience","education","platform","piety",
                 "party","gender","respondent_id", "selected", "rally")
names(cols2) = c("rally", "vote", "religion","clientelism","security","experience","education","platform","piety",
                 "party","gender","respondent_id", "selected", "rally")
names(cols3) = c("rally", "vote", "religion","clientelism","security","experience","education","platform","piety",
                 "party","gender","respondent_id", "selected", "rally")
names(cols4) = c("rally", "vote", "religion","clientelism","security","experience","education","platform","piety",
                 "party","gender","respondent_id", "selected", "rally")
names(cols5) = c("rally", "vote", "religion","clientelism","security","experience","education","platform","piety",
                 "party","gender","respondent_id", "selected", "rally")
names(cols6) = c("rally", "vote", "religion","clientelism","security","experience","education","platform","piety",
                 "party","gender","respondent_id", "selected", "rally")
names(cols7) = c("rally", "vote", "religion","clientelism","security","experience","education","platform","piety",
                 "party","gender","respondent_id", "selected", "rally")
names(cols8) = c("rally", "vote", "religion","clientelism","security","experience","education","platform","piety",
                 "party","gender","respondent_id", "selected", "rally")

bkgd <- survey[,c(grep("respondent_sex", names(survey)), grep("partisan", names(survey)), grep("respondent_sect", names(survey)),
                  grep("religious", names(survey)),grep("timer", names(survey)), grep("enumerator", names(survey)),
                  grep("low_ses", names(survey)), grep("high_ses", names(survey)), grep("years_old", names(survey)),
                  grep("household_net_income", names(survey)))]

conjoint = cbind(rbind(cols1,cols2, cols3, cols4, cols5, cols6, cols7, cols8),rbind(bkgd))

conjoint$profile = rep(c(1,2), nrow(conjoint)/2)
conjoint$round = rep(1:4,each=nrow(conjoint)/4)

#Gender
conjoint$gender <- as.factor(conjoint$gender)
conjoint$gender <- recode(conjoint$gender, "أنثى"="Woman")
conjoint$gender <- recode(conjoint$gender, "ذكر"="Man")
conjoint$gender <- relevel(conjoint$gender, ref="Man")

#Experience
conjoint$experience <- as.factor(conjoint$experience)
conjoint$experience <- recode(conjoint$experience, "لا خبرة"="None")
conjoint$experience <- recode(conjoint$experience, "5 سنوات"="5 years")
conjoint$experience <- recode(conjoint$experience, "15 سنة" = "15 years")
conjoint$experience <- relevel(conjoint$experience, ref="None")

#Sect
conjoint$religion <- as.factor(conjoint$religion)
conjoint$religion <- recode(conjoint$religion, "سنّيّ/سنّيّة"="Sunni")
conjoint$religion <- recode(conjoint$religion, "شيعيّ/شيعيّة" ="Shia")
conjoint$religion <- recode(conjoint$religion, "مسيحيّ/مسيحيّة"= "Christian")
conjoint$religion <- relevel(conjoint$religion, ref="Christian")

#Security 
conjoint$security <- as.factor(conjoint$security)
conjoint$security <- recode(conjoint$security, "في أحد الخطابات الانتخابيّة يقول/تقول المرشّح/المرشّحة: \"انّ النفوذ الايرانيّة تهدّدنا. يجب أن ندافع عن أنفسنا وشرفنا\"." = "Iranian threat")
conjoint$security <- recode(conjoint$security, "في أحد الخطابات الانتخابيّة يقول/تقول المرشّح/المرشّحة: \"انّ النفوذ الوهّابيّة تهدّدنا. يجب أن ندافع عن أنفسنا وشرفنا\"." = "Wahhabi threat")
conjoint$security <- recode(conjoint$security, "في احد الخطابات الانتخابيّة، يقول/تقول المرشّح/المرشّحة: \"المتطرّفون مناهضون للمسيحيّة. ونحن سوف نقف للدفاع عن قيمنا والحرّيّة والعدالة \"." = "Threat to Christians")
conjoint$security <- recode(conjoint$security, "في احد الخطابات الانتخابيّة، يقول/تقول المرشّح/المرشّحة: \"سأضمن أن الأمن الوطني للبنان هو الأولوية القصوى\"." = "National security")
conjoint$security <- relevel(conjoint$security, ref="National security")

#Piety 
conjoint$piety <- as.factor(conjoint$piety)
conjoint$piety <- recode(conjoint$piety, "نادرا ما يحضر/تحضر المراسم الدينيّة" = "Rarely attends religious services")
conjoint$piety <- recode(conjoint$piety, "غالبا ما يحضر/تحضر المراسم الدينيّة" = "Frequently attends religious services")
conjoint$piety <- relevel(conjoint$piety, ref="Rarely attends religious services")

#Clientelism 
conjoint$clientelism <- as.factor(conjoint$clientelism)
conjoint$clientelism <- recode(conjoint$clientelism, "يعد/تعد بالعمل  الدؤوب للناس المنتمين للدائرة الانتخابيّة" = "Promises to work hard")
conjoint$clientelism <- recode(conjoint$clientelism, "يوزّع/توزّع السلال الغذائيّة والمال خلال الانتخابات" = "Distributes food and cash")
conjoint$clientelism <- recode(conjoint$clientelism, "يؤمّن/تؤمّن العلاج الطبّيّ للناخبين"= "Arranges medical treatment")
conjoint$clientelism <- recode(conjoint$clientelism, "يساعد/تساعد الناخبين في تأمين فرص عمل لأفراد عائلاتهم" = "Helps family members get jobs")
conjoint$clientelism <- recode(conjoint$clientelism, "يساعد/تساعد الناخبين لالغاء مخالفات مواقف السيّارات" = "Assists with parking ticket cancellations")
conjoint$clientelism <- relevel(conjoint$clientelism, ref="Promises to work hard")

#Education
conjoint$education <- as.factor(conjoint$education)
conjoint$education <- recode(conjoint$education, "شهادة الثانويّة العامّة" = "High school")
conjoint$education <- recode(conjoint$education, "شهادة جامعيّة" = "University")
conjoint$education <- recode(conjoint$education,"شهادة  دراسات عليا" = "Post-graduate")
conjoint$education <- relevel(conjoint$education, ref="High school")

#Party
conjoint$party <- as.factor(conjoint$party)
conjoint$party <- recode(conjoint$party, "مرشّح/مرشّحة على لائحة الحزب المفضّل لديك" = "Favored party")
conjoint$party <- recode(conjoint$party,"مرشّح/مرشّحة على لائحة حزب يحظى بمعارضتك" = "Opposed party")
conjoint$party <- recode(conjoint$party, "مرشّح/مرشّحة على لائحة حزب لا يحظى بدعمك ولا بمعارضتك" = "Neither supported nor opposed party")
conjoint$party <- relevel(conjoint$party, ref= "Neither supported nor opposed party")

#Platform 
conjoint$platform <- as.factor(conjoint$platform)
conjoint$platform <- recode(conjoint$platform, "يعرض/تعرض خطّة مفصّلة لتعزيز خلق فرص عمل في البلاد" = "Has detailed plan for job creation")
conjoint$platform <- recode(conjoint$platform, "يعرض/تعرض خطّة مفصّلة لمعالجة مشاكل ادارة النفايات والصرف الصحّيّ في البلاد" = "Has detailed plan for waste management")
conjoint$platform <- recode(conjoint$platform, "يناقش/تناقش مشكلة ادارة النفايات والصرف الصحّيّ في البلاد" = "Discusses waste management")
conjoint$platform <- recode(conjoint$platform, "يناقش/تناقش مشكلة البطالة في البلاد" = "Discusses unemployment")
conjoint$platform <- relevel(conjoint$platform, ref="Discusses waste management")

conjoint$coethnic <- ifelse(conjoint$religion==conjoint$respondent_sect, 1, 0) 
conjoint$coethnic <- as.factor(conjoint$coethnic)

conjoint$household_ses <- as.factor(ifelse(conjoint$low_ses == 1, "Lower SES", "Higher SES"))
conjoint$income <- as.factor(conjoint$household_net_income)
levels(conjoint$income) <- c("Low Income", "Low Income",
                             "Middle Income", "Middle Income", "Middle Income", "Middle Income",
                             "High Income", "High Income", "High Income")

conjoint$ses_x_party <- NA
conjoint$ses_x_party[conjoint$household_ses == "Lower SES" & conjoint$party == "Favored party"] <- "Lower SES*Favored Party"
conjoint$ses_x_party[conjoint$household_ses == "Lower SES" & conjoint$party != "Favored party"] <- "Lower SES*Other"
conjoint$ses_x_party[conjoint$household_ses == "Higher SES" & conjoint$party == "Favored party"] <- "Higher SES*Favored Party"
conjoint$ses_x_party[conjoint$household_ses == "Higher SES" & conjoint$party != "Favored party"] <- "Higher SES*Other"
conjoint$ses_x_party <- as.factor(conjoint$ses_x_party)

#==================== conjoint analysis - main paper

# Figure 1
amce_model <- cj(conjoint,
                 selected ~ experience + gender + religion + education + party + platform + clientelism + piety + security,
                 estimate="amce", id= ~respondent_id, by= ~household_ses)

amce_model$labels <- paste(round(amce_model$estimate, digits=2),
                             paste("(",round(amce_model$std.error, digits=3),")", sep=""))
amce_model$amce_model[amce_model$labels == "0 (NA)"] <- " "

g <- plot(amce_model[amce_model$feature == "clientelism" & amce_model$level != "Promises to work hard",],
          vline=0,theme = ggplot2::theme_bw(), size=1.5) + aes(shape=household_ses, colour=household_ses)
g <- g + scale_color_manual(values=c("black","black"), na.translate = F)
g <- g + scale_shape_manual(values=c(17,19), na.translate = F)
g <- g + guides(colour = guide_legend(override.aes = list(shape = c(17,19))))
g <- g + theme(text = element_text(size=10), legend.title=element_blank(),
               legend.box.background = element_rect(colour = "black"))
g <- g + theme(axis.title = element_text(size = 9)) + theme(plot.title = element_text(hjust = 0.5))
g <- g + geom_text(size=2.2, aes(label=labels), vjust=-1.2, position = position_dodge(width=0.75), show.legend = F)
g <- g + geom_vline(xintercept=0, colour="white", linetype="longdash")
print(g)

# Figure 2
amce_model <- cj(conjoint[conjoint$household_ses == "Lower SES",],
                 selected ~ experience + gender + religion + education + platform + clientelism + piety + security,
                 estimate="amce", id= ~respondent_id, by= ~party)

amce_model$labels <- paste(round(amce_model$estimate, digits=2),
                           paste("(",round(amce_model$std.error, digits=3),")", sep=""))
amce_model$amce_model[amce_model$labels == "0 (NA)"] <- " "

g <- plot(amce_model[amce_model$feature == "clientelism" & amce_model$level != "Promises to work hard",],
          vline=0,theme = ggplot2::theme_bw(), size=1.5) + aes(shape=party, colour=party)
g <- g + scale_color_manual(values=c("black","black","black"), na.translate = F)
g <- g + scale_shape_manual(values=c(1,17,2), na.translate = F)
g <- g + guides(colour = guide_legend(override.aes = list(shape = c(1,17,2))))
g <- g + theme(legend.text = element_text(size=4), legend.title=element_blank(),
               legend.box.background = element_rect(colour = "black"))
g <- g + theme(axis.title = element_text(size = 9)) + theme(plot.title = element_text(hjust = 0.5))
g <- g + geom_text(size=2.2, aes(label=labels), vjust=-1.2, position = position_dodge(width=0.75), show.legend = F)
g <- g + geom_vline(xintercept=0, colour="white", linetype="longdash")
print(g)

# Figure 3
amce_model <- cj(conjoint[conjoint$household_ses == "Higher SES",],
                 selected ~ experience + gender + religion + education + platform + clientelism + piety + security,
                 estimate="amce", id= ~respondent_id, by= ~party)

amce_model$labels <- paste(round(amce_model$estimate, digits=2),
                           paste("(",round(amce_model$std.error, digits=3),")", sep=""))
amce_model$amce_model[amce_model$labels == "0 (NA)"] <- " "

g <- plot(amce_model[amce_model$feature == "clientelism" & amce_model$level != "Promises to work hard",],
          vline=0,theme = ggplot2::theme_bw(), size=1.5) + aes(shape=party, colour=party)
g <- g + scale_color_manual(values=c("black","black","black"), na.translate = F)
g <- g + scale_shape_manual(values=c(1,17,2), na.translate = F)
g <- g + guides(colour = guide_legend(override.aes = list(shape = c(1,17,2))))
g <- g + theme(legend.text = element_text(size=4), legend.title=element_blank(),
               legend.box.background = element_rect(colour = "black"))
g <- g + theme(axis.title = element_text(size = 9)) + theme(plot.title = element_text(hjust = 0.5))
g <- g + geom_text(size=2.2, aes(label=labels), vjust=-1.2, position = position_dodge(width=0.75), show.legend = F)
g <- g + geom_vline(xintercept=0, colour="white", linetype="longdash")
print(g)

#==================== appendix

# Figure A.2
amce_all <- cj(conjoint,
               selected ~ experience + gender + religion + coethnic + education + party + platform +
                 clientelism + piety + security, estimate="amce", id= ~respondent_id)

amce_all$labels <- paste(round(amce_all$estimate, digits=2), paste("(",round(amce_all$std.error, digits=2),")", sep=""))
amce_all$labels[amce_all$labels == "0 (NA)"] <- " "

plot(amce_all,
     vline = 0, theme=ggplot2::theme_bw(), size = 1.5, legend_title = " ") +
  theme(text = element_text(size=10)) +
  facet_grid(feature~., scales="free") +
  scale_color_manual(values=c("black","black","black","black","black","black","black","black","black","black")) + 
  geom_text(size=1.4, aes(label=labels, vjust=-1.2, hjust=0.4)) +
  guides(color = FALSE)

# Figure A.3
amce_model <- cj(conjoint,
                 selected ~ experience + gender + religion + education + party + platform + clientelism + piety + security,
                 estimate="amce", id= ~respondent_id, by= ~income)

amce_model$labels <- paste(round(amce_model$estimate, digits=2),
                           paste("(",round(amce_model$std.error, digits=3),")", sep=""))
amce_model$amce_model[amce_model$labels == "0 (NA)"] <- " "

g <- plot(amce_model[amce_model$feature == "clientelism" & amce_model$level != "Promises to work hard",],
          vline=0,theme = ggplot2::theme_bw(), size=1.5) + aes(shape=income, colour=income)
g <- g + scale_color_manual(values=c("black","black","black"), na.translate = F)
g <- g + scale_shape_manual(values=c(15,17,19), na.translate = F)
g <- g + guides(colour = guide_legend(override.aes = list(shape = c(15,17,19))))
g <- g + theme(text = element_text(size=10), legend.title=element_blank(),
               legend.box.background = element_rect(colour = "black"))
g <- g + theme(axis.title = element_text(size = 9)) + theme(plot.title = element_text(hjust = 0.5))
g <- g + geom_text(size=2.2, aes(label=labels), vjust=-1.2, position = position_dodge(width=0.75), show.legend = F)
g <- g + geom_vline(xintercept=0, colour="white", linetype="longdash")
print(g)

# Figure A.4
diff_amce <- cj(conjoint,
                selected ~ experience + gender + religion + education + platform + clientelism + piety + security,
                estimate="amce_diff", id= ~respondent_id, by= ~household_ses)
diff_amce$labels <- paste(round(diff_amce$estimate, digits=2), paste("(",round(diff_amce$std.error, digits=2),")", sep=""))

g <- plot(diff_amce[diff_amce$feature == "clientelism",],
          vline=0,theme = ggplot2::theme_bw(), size=1.5) + aes(shape=household_ses, colour=household_ses)
g <- g + scale_color_manual(values=c("black"), na.translate = F)
g <- g + scale_shape_manual(values=c(19), na.translate = F)
g <- g + guides(colour = guide_legend(override.aes = list(shape = c(19))))
g <- g + theme(text = element_text(size=10), legend.position = "none") +
  xlab("Estimated AMCE Difference (Low SES - High SES)")
g <- g + theme(axis.title = element_text(size = 9)) + theme(plot.title = element_text(hjust = 0.5))
g <- g + geom_text(size=2.2, aes(label=labels), vjust=-1.2, position = position_dodge(width=0.75), show.legend = F)
g <- g + geom_vline(xintercept=0, colour="white", linetype="longdash")
print(g)

# Figure A.5
conjoint$party <- relevel(conjoint$party, ref="Opposed party")

diff_amce <- cj(conjoint[conjoint$household_ses == "Lower SES",],
                selected ~ experience + gender + religion + education + platform + clientelism + piety + security,
                estimate="amce_diff", id= ~respondent_id, by= ~party)
diff_amce$labels <- paste(round(diff_amce$estimate, digits=2), paste("(",round(diff_amce$std.error, digits=2),")", sep=""))

g <- plot(diff_amce[diff_amce$feature == "clientelism" &
          diff_amce$BY != "Neither supported nor opposed party - Opposed party",],
          vline=0,theme = ggplot2::theme_bw(), size=1.5) + aes(shape=party, colour=party)
g <- g + scale_color_manual(values=c("black"), na.translate = F)
g <- g + scale_shape_manual(values=c(19), na.translate = F)
g <- g + guides(colour = guide_legend(override.aes = list(shape = c(19))))
g <- g + theme(text = element_text(size=10), legend.position = "none") +
  xlab("Estimated AMCE Difference for Low SES (Favored Party - Opposed Party)")
g <- g + theme(axis.title = element_text(size = 9)) + theme(plot.title = element_text(hjust = 0.5))
g <- g + geom_text(size=2.2, aes(label=labels), vjust=-1.2, position = position_dodge(width=0.75), show.legend = F)
g <- g + geom_vline(xintercept=0, colour="white", linetype="longdash")
print(g)

# Figure A.6
diff_amce <- cj(conjoint[conjoint$household_ses == "Higher SES",],
                selected ~ experience + gender + religion + education + platform + clientelism + piety + security,
                estimate="amce_diff", id= ~respondent_id, by= ~party)
diff_amce$labels <- paste(round(diff_amce$estimate, digits=2), paste("(",round(diff_amce$std.error, digits=2),")", sep=""))

g <- plot(diff_amce[diff_amce$feature == "clientelism" &
                      diff_amce$BY != "Neither supported nor opposed party - Opposed party",],
          vline=0,theme = ggplot2::theme_bw(), size=1.5) + aes(shape=party, colour=party)
g <- g + scale_color_manual(values=c("black"), na.translate = F)
g <- g + scale_shape_manual(values=c(19), na.translate = F)
g <- g + guides(colour = guide_legend(override.aes = list(shape = c(19))))
g <- g + theme(text = element_text(size=10), legend.position = "none") +
  xlab("Estimated AMCE Difference for Low SES (Favored Party - Opposed Party)")
g <- g + theme(axis.title = element_text(size = 9)) + theme(plot.title = element_text(hjust = 0.5))
g <- g + geom_text(size=2.2, aes(label=labels), vjust=-1.2, position = position_dodge(width=0.75), show.legend = F)
g <- g + geom_vline(xintercept=0, colour="white", linetype="longdash")
print(g)
